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ABSTRACT 

We have performed three-dimensional magneto-hydrodynamical simulations 
of stellar accretion disks, using the PLUTO code, and studied the accretion of 
gas onto a Jupiter-mass planet and the structure of the circumplanetary gas 
flow after opening a gap in the disk. We compare our results with simulations 
of laminar, yet viscous disks with different levels of an a-type viscosity. In all 
cases, we find that the accretion flow across the surface of the Hill sphere of 
the planet is not spherically or azimuthally symmetric, and is predominantly 
restricted to the mid-plane region of the disk. Even in the turbulent case, we 
find no significant vertical flow of mass into the Hill sphere. The outer parts of 
the circumplanetary disk are shown to rotate significantly below Keplerian speed, 
independent of viscosity, while the circumplanetary disk density (therefore the 
angular momentum) increases with viscosity. For a simulation of a magnetized 
turbulent disk, where the global averaged alpha stress is olmhd = 10~ 3 , we find 
the accretion rate onto the planet to be M ~ 2 x 10~ 6 Mjyr~ 1 for a gap surface 
density of YlgcmT 2 . This is about a third of the accretion rate obtained in a 
laminar viscous simulation with equivalent a parameter. 

Subject headings: planet accretion, magneto-hydrodynamics, protoplanetary disks 
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INTRODUCTION 



Gas giants can form as a result of core formation by planetesimal accretion followed by 



formation of the gas envelope through accretion from the circumstellar disk (IMizuno 



Pollack et al. 



1980h . 



(119961 ) distinguished three phases of this process, using numerical simulations 



of core accretion and envelope evolution. A first phase marked by the fast accretion of solids 



19821 ): a 



onto a core until the feeding zone of the planet is mostly evacuated ( jStevensonl 
second phase where gas and solid accretion is low and constant; a third stage when the core 
mass equals the envelo pe mass, leadi ng to the contraction of the envelope and the onset of 



runaway gas accretion (IMizuno 



19801 ). Migration of the planet might allow for an extension 



of the feeding zone, while gap formation might lead to a mass limit for gap opening planets 



( lAlibert et al.l 120051 ). In the outer parts of the disk, giant planets could form as a result 



of the collapse of a gravitationally unstable disk clump. This mechanism requir es a yery 



massive disk that can cool effect i vely on timescales o f a few local orbital periods (IBoss 



1997 



Mayer et al. 



2002 



Rafikov 



2005 



Durisen et al. 



2007!) 



Monte Carlo models of planet formation produce synthetic populations of extrasolar 
planets, with a large diversity of initial conditions. These models have bee n successful 



i n rep r oducing key features of t he observed dis t ribut ion of exoplanets (e.g. 



(120081 ): 



Mordasini et al. 



(12009b 



Benz et al. 



aj); llda and Linl ( 120081 )). The calculations usually include 



one- dimensional disk evolution and core/envelope structure models. The accretion of 
planetesimals and gas onto an already formed proto-core is included, using a given 



presc ription for the accretion rates of gas and rocky materials onto the planet (lAlibert et al. 
2004 ). For this reason, an accurate estimation and parameterization of the accretion rates of 
gas onto planets for a variety of conditions is necessary to correctly calculate the formation 
time and the limiting mass of the giant planetary population. 



1 Currently at University of Chicago 
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So far, the accretion of gas onto planets has been modeled using two different approaches. 
On one side, one dimensional models have been used to calculate the radial structure of 
the envelope and the accretion onto a rocky core. These models can include effects such as 



the dust opacity of the envelope, the release of e nergy of infa. 



mg_p 



2000 



anetesimals into the 



Hubickyj et al. 



2005|). 



envelope and the thermal feedback of the planet (llkoma et al. 
These models can include the disk evolution only in a restricted way, and assume a certain 
model of the outer envelope (as a boundary condition) that is spherically symmetric. 

On the other hand, two/three-dimensional simulations of accretion disk with accreting 
planets aim to estimate the structure of the flow around the planet and how much mass 
the disk is capable of feeding to the planet. However, most of these simulations do not 



include the radiative feedback 



Tanigawa and Watanabe 



rom t he planet and a detailed model of the inner envelope. 



J2002J) and 



D'Angelo et al. 



2003aU bl) used high-resolution 



two-dimensional simulations to study the detailed flow pattern around the planet, and the 
circumplanetary disk. They showed that inside the planet's Roche lobe, accretion in the 
circumplanetary disk is powered mainly by energy dissipation of circulating matter at the 
spiral shock. Outside the Roche lobe, gas flows onto the planet through "accretion bands" 
located between the horseshoe flow and the passing-by flow (although the detailed structure 
depends strongly on the sound speed). The accretion timescale 



Mp 



(1) 



has been measured to be around 10 4 — 10 5 yr, and the accretion rate of a Jupite r-mass planet 



has been found to 



1999 



je of the order of 10 5 Mj/yr on a disk w ith = O.O1M (IBryden et al. 



Lubow et al 



1999 



Kley et al 



2001 



Bate et al 



20031 ) . Three-dimensional simulations 



including radiation transfer have found simila r accretion rates and h ave shown the formation 



of a thick (H/r « 0.5) circumplanetary disk (IKlahr and Kley 



20061 ) 



The properties of the circumplanetary disk have been the subject of various numerical and 
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semi-analytical studies. lAyliffe and Bate! ( 120091 ) used radiation-hydrodynamics numerical 
simulations to study the planet-disk system, and derived important properties such as 
the power law dependence of density with radius, and the scale height of the stellar disk. 
Their calculations show that the disk density falls with radius as r~ 2 to r~ 3 / 2 , and the 



Ward and Canup 



f hoiol ) 



disk is fairly thick, with H/r > 0.2, with no significant flaring, 
developed a model for the circumplanetary disk that includes the contraction of the planet 
envelope and the properties of the inflowing gas material from the disk. They find various 
phases of gas accretion and collapse where the circumplanetary disk can either provide 
mass to the planet or act as a "spin-out disk" where the mass flux is outwards. The inner 
circumplanetary disk rotates Keplerian up to a small fraction of the Hill sphere and then 
rotation becomes sub-Keplerian in the outer region. The issue of the circumplanetary disk 
size, or the truncation mechanism of the circumplanetary disk, is also crucial for the fate 
of the solid material. The inner Keplerian rotating thick disk is truncated at a fraction of 
the Hill sphere (oc 0.4 r/,,), but this seems to depend on the viscou s mechanism acting in the 



circumplanetary disk (IMartin and Lubow 



2011 



Bate et al. 



2003|). 



Recent work by iTanigawa and Ikomal (j2007l ) suggests a different picture for the accretion 
flow in an inviscid isothermal circumplanetary disk. In this case, the bulk of the accretion 
does not occur in the disk mid-plane region, but instead at high latitudes in the planet 
Hill sphere. In the equatorial region, an outflow is present thr ough the Lagran g ian p oints 



Klahr and Klevl fl2006f ) in 



LI and L2. This flow structure is similar to that observed by 
radiative non-viscous simulations. 

In this paper we perform numerical simulations to study the accretion rate of gas onto giant 
gap-opening planets and the circumplanetary disk in turbulent magnetized disks using, for 



the first time, three-dimensional global stellar disk simulat ions. The turbulence in t 



re disk 



is generated by the magneto-rotational instability (MRI) (jBalbus and Hawley 



19911). A 



minimum degree of ionization leading to a good coupling between the gas and the magnetic 
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field in the stellar disk is required for the MRI to function. In the cold and dusty mid-plane 
regions the ionization will be low, possibly producing a dead-zone, while on the outer 
parts of the disk, cosmic rays can penetrate allowing; for the n ecessary ionization degree 



( iTurner et al. 



2010 



Dzyurkevich et al. 



2010 



Flock et al. 



201 lh . A planet gap, however, 



might lower the local density to a sufficient level to allow for ionization of the mid-plane 
regions, leading to an MRI active zone. Therefore, this type of turbulence could conceivable 
be still present in the planet corotation region. 

We compare the accretion rate onto a planet in an MRI-turbulent disk with that in a 
viscous laminar disk and we examine the accretion structure and mass inflow into the Hill 
sphere, where material passes through the circumplanetary disk to be accreted by the 
planet. The viscosity provided by the turbulence in the stellar disk (either modeled by an 
a prescription or by "real" magnetic turbulence) is acting indirectly in the circumplanetary 
disk. There is no extra viscosity applied to the circumplanetary disk, nor is this disk MRI 
unstable on its own. All the effects we observe are due to the viscosity of the circumstellar 
disk, be this a laminar kinematic viscosity, or a magnetic turbulent transport. Note that 
the MRI is not resolved in the circumplanetary disk, and the question of what mechanism 
drives accretion in this disk is still unresolved. The effective transport comes from the 
transport of angular momentum in the stellar disk. We also concentrate strictly on the 
phase when a gap has been cleared in the disk, and study the efficiency of the stellar disk 
in providing mass to the forming protoplanet. The choice of simulating a locally isothermal 
disk is justified for simulations with a Jupiter mass planet since at this planet mass the 
accretion rate approach es the isothermal result and is independent of the opacity treatment 



(Ayliffe and Bate 



2009). Additionally, for this mass range, the accretion rate onto the 



planet depends on the ability of the circumstellar disk to provide material, which is the 
effect we study in this paper. 



The paper is organized as follows. In Section 2, we describe the computational setup, 
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boundary and initial conditions, and the parameters we use in our simulations. We also 
describe the prescription for the mass accretion onto the planet. In Section 3.1, we present 
our results on the three-dimensional structure of the accretion flow into the Hill sphere. 
Section 3.2 contains the results on the mass accretion rates for the different simulations, 
while Section 3.3 shows additional details of the accretion flow around the planet. Finally, 
we discuss and summarize our results in Sections 4 and 5. 



2. COMPUTATIONAL SETUP 



Simulations were per formed using the finite volume fluid dynamics code PLUTO 



( jMignone et al. 



20071 ). In the code, time stepping is done using a second order Runge Kutta 
scheme, while the spatial integration is performed using linear interpolation through the 
second order TVD scheme. The Riemann fluxes are computed using the HLLC and HLLD 



solvers for the HD and MHD cases, respectively. The co de uses the Constrained 



method for preserving a divergence-free magnetic field (I Gardiner and Stone 



numerical setup for the MHD case follows the setup presented in 



Flock et al. 



ransport 



2005). The 



(120101 ). 



We use spherical coordinates (r, 6, <p) and the domain is given by r G [1, 10] (in units of AU), 
6 G [vr/2-0.3, tt/2 + 0.3] and <p G [0, 2tt]. The grid resolution is (N r , N g , = (256, 128, 256) 
and it is centered in the center of mass of the planet-star system. 

The circumstellar gas disk is initially in sub-Keplerian rotation around a solar mass star. 
The azimuthal velocity is given by 



v</> = \fk -c 2 s (a-2b), 



(2) 



where Vk is the Keplerian velocity and a = 3/2 and b = 0.5 are the exponents of the radial 



power law distribution of the density p oc r a and sound speed c s = Co(r sin 6) . The initial 



8 



density distribution is given by 



p(r,6) = (rsintf)- 3 / 2 



cxp 



( 



sin 9 — 1 




(3) 



The disk is described by a locally isothermal equation of state P = c 2 s p. The ratio of the 
pressure scale height h to the radial coordinate of the disk h = c s /v k is taken to be a 
constant such that h = H/{r sin 0) = 0.07. 

The gravitational potential of the planet is given by a softened point-mass potential 



where e is the softening parameter, needed to avoid numerical divergence near the 
position of the planet. For all the simulations e is set to be a fraction of the Hill radius 
e = kr p (M p /3) 1 ^ with k = 0.3. Distances are given in units of r = 1AU, density is given 
in units of p = 1 x 10~ 12 gcm~ 3 , and velocity is given in units of Keplerian speed at 1AU, 
v = Vk(lAU). The density has been scaled such that the total disk mass is 0.01M stor . 
Magnetic fields are given in units of B = ^4ttp Vq. The equations of motion of the planet 
are solved at each timestep with a leap frog integrator. 



The boundary conditions for the velocities and magnetic field are periodic in the vertical 
(9 boundary) and azimuthal directions and reflective in the radial direction, except for the 
transverse magnetic field component, which reverses its sign at the radial boundary. Buffer 
zones are defined at the radial boundaries to avoid boundary effects, where for 1 < r < 2 
the magnetic resistivity is given by i] = 2 x 10~ 4 (2 — r) and for 9 < r < 10 the resistivity is 




GM 



(4) 




2.1. Boundary conditions 



77 = 1 x 10- 4 



(r-9). 
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2.2. Initial conditions, gap opening and viscosity 

Since our aim is to find the accretion rate onto the planet only in the stage when the 
gap is cleared, the planet is allowed to accrete gas after 100 orbital periods at 5AU 
have elapsed (see Figured]). At this stage, a gap has been cleared, and the density has 
been reduced by more than 95%. For the hydrodynamical simulations, viscosity is added 
explicitly as a source term in the momentum equation. We use an a-type kinematic 
viscosity in the circumstellar disk which is given by v = ac s H, where a takes values of 
2 x 10" 3 , 1 x lO" 3 , 5 x lO" 4 , and 1 x 10" 4 . 

2.3. Accretion prescription 

The accretion onto the planet is modeled by removing a fraction of the mass inside a 
portion of the Hill sphere at each time step. At each timestep the new density p is given by 

P(r)= (l-^)p(r). (5) 

The accreted mass by the planet in timestep At is AM = [p{r)Att' a 1 )r 2 sin{9)drd9d(f). 
The planet accretion rate for timestep At is calculated as the accreted mass divided by 
the timestep AM/ At. The factor t a represents the accretion timescale in which the Hill 
sphere is emptied if there was no gas flowing in from the disk. This is chosen to be 
t a = lfi^^ inside the inner half of the Hill sphere and t a = 2^^^ in the outer half of the 
Hill sphere (where fi^y is the Keplerian angular frequency at 1AU). A density floor is 
applied to the simulations with magnetic fields, where the density is not allowed to drop 
below 10 -19 gcm~ 3 . Nevertheless, the density in the simulation never reaches this value. 
The magnetic field is not modified as the density is reduced inside the Hill sphere in order 
to preserve a divergence-free field. The accretion rate has been shown to be dependent 
on the accretion radius (the distance from the planet up to which mass is removed) and 
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Fig. 1. — Initial conditions before the planet starts accreting for the laminar viscous disk 
simulations and the MHD simulation. The gap in the MHD simulation is found to be wider 
as compared to all the viscous simulations. The figure shows the density after 100 orbital 
periods have elapsed. 



on the accretion timescale parameter t a . 



Tanigawa and Watanabd (120021 ) showed that the 



accretion radius should be small (~ O.lr^) and the accretion timescale should be on the 
order of the orbital period, in order to obtain converged results. Because of our lower 
resolution, we take most of the mass from within the inner half of the Hill sphere (our 
resolution element is 0.035r or 0.1 r^j ). This prescript ion has also been used in previous 



studies of gas accretion and migration by giant planets (IKley et al. 



20011 ). We also checked 
that our results remain valid if one restricts the accretion radius down to OAr^. We have 
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also verified that we obtain the same results if we extend the accretion timescale to t a = 10. 



3. RESULTS 

3.1. Structure of the envelope and mass inflow 

In this section we study the structure of the density and the inflow of mass in/into 
the Hill sphere. The density structure of the disk around the planet can be seen in 
Figure [2j The color and contour lines represent the density. The circumplanetary disk 
is shown for the viscous simulations with a = 2 x 10~ 3 , 1 x 10~ 3 , 5 x 10~ 4 , and for the 
magnetized disk simulation (bottom right panel). The horizontal density contours signal 
a disk with a scale height increasing linearly with radius from the planet, and with no 
flaring. Additionally, the density drops with radius away fro m the planet app r oxima tely as 



Avliffe and Batel (120091 ) using 



~ r~ 3 / 2 . This is consistent with the density profile found by 
higher-resolution radiation-hydro simulations. The find the disk profile to be described by 
a power law between r~ 2 and r _3//2 . Figure [3] shows the azimuthal velocity and density 
of the circumplanetary gas. The velocity profile is consistent with rotation around the 
planet, and the density profile increases inwards in all cases. We note that the rotation of 
the circumplanetary gas is significantly sub-Keplerian (with respect to Keplerian motion 
around the planet, not the star). The azimuthal velocity of the gas in the disk of the planet 
rotates at less than 50% Keplerian speed. The velocity structure is also not significantly 
sensitive to the viscosity, while the density (and therefore the angular momentum) increases 
with viscosity. The higher density in the large viscosity runs might be due to the higher 
mass transfer rate from the circumstellar disk to the circumplanetary disk. The larger the 
viscosity, the larger the amount of mass that the stellar disk can feed to the planet disk, 
therefore increasing the density surrounding the planet, and possibly the accretion rate of 
gas onto the planet. 
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The distribution of angular momentum in the circumplanetary disk is shown in Figure HJ 
The left plot of this figure shows the angular momentum scaled with the initial density, 
while the right plot shows the specific angular momentum. Both quantities are calculated 
with respect to the planet. The angular momentum distribution increases in the inner 
part of the disk, and levels off at ~ 0.2 — 0.3rh- This is c onsist ent with the definition of 



the edge of the planet disk proposed by 



Ayliffe and Bate! (120091 ). since in a Keplerian disk 



angular momentum increases monotonically outwards (however, at this distances we start 
to reach the limit of our grid resolution which is O.lHillradii in the radial direction). 
This criterion could define the edge of the Keplerian rotating inner disk, before rotation 
becomes sub-Keplerian in the outer parts of the disk. However, we also note that we do not 



obtain Keplerian rotation in t 



limit. 



Tanigawa et al. 



lis inner part, although this might be affected by a resolution 



Avliffe and Batd (120091 ) obtained a value of ~ 0.3r h for the cutoff radius, and 



(120121 ) obtained similar values. This value is also similar to the truncation 



radius defined by the radius where tidal effect s by the central star beco me important enough 



to disturb the inner Keplerian rotating disk (jMartin and Lubow 



201lh . 



We calculate the mass flux through the surface of the Hill sphere. The mass flux is given 



by pVi 



inflow 



pv • (VF/|VF|), where F(r) 



is the equation describing 



the surface of the Hill sphere. The mass flux is plotted in Figures [5] and [6] for four runs. 
These quantities have been averaged in time. We plot the surface of the Hill sphere using 
an ellipsoidal projection where the three-dimensional structure can be observed. The center 
of the ellipse corresponds to the point in the Hill sphere that is most distant from the star. 

The density and mass flux are much larger in the case with the higher viscosity 
(a = 2 x 10~ 3 ), and are the lowest in the turbulent magnetized case. In all cases, all the 
net flux through the surface is inflowing, and we see no significant amount of gas entering 
the Hill sphere from the vertical direction above and below the mid-plane. The bulk of the 
accretion at this distance from the planet is confined to the mid-plane "region". Figure [7] 
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shows the vertical structure of the mass flux averaged over the azimuthal direction. 




0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 

(r-r„)/r s (r-r p )/r h 



Fig. 2. — Density structure of the disk around the planet (in the radial and vertical direction). 
The planet is located at the origin of the coordinate system. Color and contour lines show 
the logarithm of the density. Top left: a = 5 x 10~ 4 . Top right: a = 10~ 3 . Bottom left: 
a = 2 x 1CT 3 . Bottom right: MHD simulation. 



3.2. Mid-plane flow structure in the vicinity of the planet 

In this section we study the dependence of the flow structure on the viscosity in the 
simulations and the effect of turbulence on the accretion flow onto the planet. Figure [8] 
shows the radial mass flux and the pressure (and magnetic pressure b 2 /(8n)) at a distance 
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Fig. 3. — Left: Azimuthal velocity of the circumplanetary gas (relative to Keplerian speed). 
The dashed line shows the Keplerian speed profile relative to the planet, that has been 
reduced by 50%. Right: Density of the circumplanetary gas for the different runs. These 
figures show radial profiles in a line joining the star and planet. 




Fig. 4. — Left: Angular momentum of the circumplanetary gas (normalized to initial densi- 
ties). Right: Specific angular momentum of the circumplanetary gas for the different runs. 
These figures show radial profiles in a line joining the star and planet. 




Fig. 5. — Mass flux through the surface of the Hill sphere for the viscous laminar runs with 
a = 5 x lCT 4 (left) and a = 1 x 10~ 3 (right). The mass flux is given in units of Mjyr~ 1 S~ 1 , 
where the quantity S is the area of the grid cell in the Hill sphere given by S = t^AOhs^-^hs- 
The center of the ellipse corresponds to the point in the Hill sphere that is most distant from 
the star and points away in the radial direction. 




Fig. 6. — Mass flux through the surface of the Hill sphere for a = 2 x 10~ 3 (left) and for 
the MHD-turbulent run(right). The mass flux is given in units of Mjyr~ l S~ 1 , where the 
quantity S is the area of the grid cell in the Hill sphere given by S = t^AOhs^hs- The 
center of the ellipse corresponds to the point in the Hill sphere that is most distant from the 
star and points away in the radial direction. 

of ±2r^ from the planet position. In the mid-plane, the thermal pressure exceeds the 
magnetic pressure by a factor of the order of 10 2 . There is radial inflow of gas coming from 
both sides of the planet (although inflow is not spherically symmetric in the mid-plane 
and is not maximum in the line joining the star and planet), and this is dominant in the 
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Fig. 7. — Vertical structure of the mass inflow pv in fi ow into the Hill sphere. The coordinate 
6rh refers to the polar angle in the frame of the planet. The quantity pv in fi ow has been 
azimuthally averaged (with respect to the Hill sphere). The quantity S is the area of the 
grid cell given by S = rlA6 RH A(j) RH . 

high viscosity case. As the viscosity gets lower, the density and mass flux decrease. The 
magnetized case shows radial mass inflow comparable to the two cases with the lower 
a viscosity, or lower. The radial profile of the thermal pressure is similar for all cases, 
although the pressure across the Hill sphere decreases with viscosity. At either side of the 
planet, following the spiral arms, there are bumps of high magnetic pressure. In this region, 
the magnetic pressure is now comparable to the thermal pressure. 



Figure |9] shows the density in the mid plane and the radial velocity for the laminar viscous 
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simulation with a = 2 x 10 3 . Overp lotted is the mid-plane vector field of the velocity. 
Figure [TU] shows the same quantities for the ma gnetized turbulent circumstella r disk 



simulation. In agreement with previous studies ( ITanigawa and Watanabe 



2002), we find 



that the gas accreting onto the planet comes from a flow between the open pass-by flow and 
the gas that is orbiting in horseshoe orbits at corotation. This comes from both sides of the 
planet, as can be seen in the velocity field in the upper right and lower left part of Figures 
[9] and [TO]. Material enters the Hill sphere through these channels, as it is also seen in the 
mass flux at the surface of the Hill sphere in Figures [5] and [6] Inside the Hill sphere, the 
gas accretes onto the planet through the circumplanetary disk. Two-dimensional inviscid 
simulations show accretion in the planet disk that is powered by a spiral shock where 
gas looses energy and spirals inwards. This is not observed in our simulations, since in 
three-dimensions shocks are smoothed out and the resolution in the disk is not sufficient. 
Outside the Hill sphere, we see the spiral arm structure that forms the bow shock (see 



Tanigawa and Watanabd (120021 ); 



D'Angelo et al. 



(j2003bl )). although the shock is smeared 



out by viscosity and turbulence in our simulations. 

In the case of the magnetized turbulent run, the velocity structure around the planet is 
much less uniform in comparison with the laminar viscous run. This is due to small scale 
turbulence and the non-uniformity of the magnetic field. However, the mean velocity 
structure, averaged in time, is similar for both cases as is shown in Figures 191 and [TU1 



3.3. Gas accretion rates and planet growth time after gap-opening 



The cumulative accreted mass and the accretion rate onto the planet are shown in Figures 
[T2l and [HI for the laminar viscous simulations and the magnetized simulation. We will refer 
to en mh d to denot e the alpha stress that is me asured in the magnetized simulation, and 
to a to denote the Ishakura and Sunyaevl ( 1973 ) viscosity parameter that is chosen for the 
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Fig. 8. — Left: Radial mass flux for the different runs. Right: Pressure for the different 
runs. The dashed line shows the magnetic pressure (multiplied by a factor of 50) for the 
magnetized case. 




-0.2 




Fig. 9. — Left:Density (in units of 10~ 12 grcm~ 3 ) in the mid-plane for the laminar viscous 
simulation with a = 2 x 10~ 3 . Right: Radial velocity (in units of v^ilAU)) in the mid- 
plane for the same simulation. The overplotted vector field shows the velocity field in the 
mid-plane. 

viscous laminar simulations. The measured alpha is calculated as sum of the contributions 
to the internal stress by the Reynolds Tr =< pSv^Sv? > and Maxwell magnetic stress 
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Fig. 10. — Left:Density (in units of 10~ 12 grcm~ 3 ) in the mid-plane for the MHD simulation. 
Right: Radial velocity (in units of Vk(lAU)) in the mid-plane for the MHD simulation. The 
overplotted vector field shows the velocity field in the mid-plane. 

Tm =< B^Br/lAirP) >, where the brackets denote the mean value over the computational 
domain (excluding magnetic buffer zones) and P is the pressure. The total stress is 
additionally normalized by the density. The velocities <5t> r ,0 = v r ^— < v r ^ > denote 
deviations from the mean value. 

We will first discuss the laminar simulations. The largest planet accretion rate is 
obtained for the viscous simulation with a = 2 x 10~ 3 , as it is expected since the 
stellar disk accretion rate is proportional to the viscosity in the more simple analytical 
approximation. The limit of the lowest viscosity that the code is able to resolve above 
numerical dissipation effects is a ~ 10~ 4 . Additionally, at this low viscosity, our simulation 
time is not able to cover the viscous evolution, since the viscous timescale is given by 
r msc = H 2 /u = H 2 /(aH 2 Q) = (orfi)- 1 . 

The magnetic case shows an interesting behavior. We find the average planet accretion 
rate in the magnetic case to decrease drastically during the first stages of the simulation. 
It seems to stabilize at a value of the accretion rate similar all the viscous simulations, but 
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still presents more variability than the purely viscous cases. In the magnetized case, the 
planet accretes gas at a rate which is around a third of the accretion rate in an disk with 
a = 1 x 10~ 3 . For the turbulent magnetized simulation, the global and time-averaged oimhd 
is equal to olmhd = 1 x 10~ 3 . The global average Maxwell stress is atMHD,Max = 2 x 10~ 4 . 
Due to the presence of the giant planet, the Reynolds stress dom inates over the Maxwell 



stress by a factor of 2 to 3 ([Winters et al. 



2003 



Uribe et al. 



201ll ). However, the effective 



viscosity provided by the turbu 



circumstellar disk (IFlock et al. 



ence in the mid-plane is less than in the upper layers of the 



20111 ). Small scale turbulent structures in the mid-plane 
might not be well resolved and the a parameter measured for magnetic turbulence measures 
large scale transport. Nevertheless, the effective viscosity in the mid-plane should be 



comp arable to the one in the viscous laminar simulation with a = 1 x 10 4 (IFlock et al. 



201lh . 



Figure [12] shows the total mass accreted by the planet starting at a time when the planet 
accretion is switched on. There is an initial rapid rise due to the material that has been 
accumulated in the Hill sphere during the previous evolution of 100 orbital periods. After 
this stage, the planet has consumed the "excess" of material, and accretes at a rate by 
which the circumplanetary disk can provide material. It can be seen in Figure [12] that 
even though the initial conditions are slightly different (some simulations have more gas 
accumulation around the planet depending on viscosity, as seen in Figured]), after the initial 
phase is passed, the planet accretion tends to relax to more or less steady state values. 



Previousl y iTanigawa and Watanabd (120021 ) and 



4xl0 4 yr. 



Kiev et al. 



Lubow et al. 



(119991 ) found growth times of 



( 1200 ll ) found gas accretion rates onto the planet of M = 6xl0 5 Mj/y r 
for Ju piter mass planets, which indicates a growth time of 2 x 10 4 yr. (lAyliffe and Bate 



2009bl ) calculated accretion rates for various values of the planet mass, using radiation- 
hydrodynamics simulations. For a planet with mass equal to the mass we use in our study 
(a Jupiter mass), they find an accretion rate onto the planet of ~ 3 x 10 _5 M/ up ?/r _1 . 
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Fig. 11. — Mass accretion rate onto the planet for the four viscous simulations and the 
magnetized simulation. The yellow line shows the MHD case. The colored dashed lines 
show the mean value of each simulation (MHD rates are averaged after 30 orbits). 

We find a mass accretion rate of M ~ 9 x 10 -7 Mj/yr for the laminar viscous case with 
a = 2 x 10~ 3 . For the magnetized simulation we find M w 2 x 10~ 7 Mj/yr. These rates are 
measured only when the planet has already cleared a gap around its orbit, contrary to the 
previous studies. Our surface density at the position of the planet is therefore very small, as 
compared to previous studies. We measure a surface density after gap opening of l.2gcm~ 2 
at the position of the plane t. Scaling our accretion rates to the surface density value used 



by 



Ayliffe and Bate! (j2009bl ). the accretion rate we measu re is 3 x 10 5 Mj/yr, which is 



similar to what was obtained by lAyliffe and Batd (l2009bl ). The scaling of the accretion 
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Fig. 12. — Cumulative mass accreted by the planet for the four viscous simulations and the 
magnetized simulation. 

rate is possible, since the quantity should vary linearly with surface density (in the simple 
viscous model). 



4. DISCUSSION AND CONCLUSIONS 

Figure [13] shows the planet accretion rate as a function of the circumstellar disk a for the 
laminar viscous runs. The rate obtained in the magnetized run is shown as a dotted line 
since a is not constant in this case, and the diamond symbol signals the global stress at 
the beginning of the simulation. For the magnetized case, after 100 orbital periods, the 
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turbulence has decayed as was seen in the simulations presented in lUribe et al.l ( 1201 ll ). 
The effective global averaged stress coming from the turbulence at the beginning of the 
simulation (right before planet accretion starts) is olmhd = 1 x 10~ 3 . The Maxwell stress 
at the beginning of the run is olmhd, Max = 2 x 10~ 4 . In the magnetic case, the measured 
planet accretion rate is smaller than in all laminar viscous runs and it is a third of the 
accretion rate in the viscous run with a ~ 1 x 10~ 3 . This can be attributed to the fact 
that turbulent transport is achieved mainly at the large scales, while the effective viscosity 
provided by the turbulence at the small scales is not represented by the global value of the 
measured olmhd- There also remains the question of how well small-scale turbulence is 
resolved in our simulations. It is possible that higher resolution allows for smaller structure 
turbulence and transport on smaller scales. 



The mass accretion into the Hill sphere happens 



each side of the planet ( ITanigawa and Watanabe 



along two channels or bands located at 



20021 ). Closer to the radial location of 



the planet, material cannot flow in, and instead performs a U-turn, since it looses its 
angular momentum rapidly as it approaches the planet. Radially away from the planet, the 
gravitational torque of the planet is not strong enough to pull the material in fast enough, 
and instead gas orbits passing the planet. The planet accretion flow lies between these two 
regions. In all cases, we find that the accretion of gas into the Hill sphere is not spherically 
symmetric, nor azimuthally symmetric. In all cases, the flow through the vertical direction 
is negligible and the flow is restrained to low latitudes. In our simulations, the Hill radius 
is approximately equal to the pressure scale height of the circumstellar disk. 



Ayliffe and Batd (120091 ). who performed 



Our results are consistent with results presented by 
hydrodynamic radiation simulations of the planet-disk system. They show that the inflow 
of the planet accretion flow is confined to the mid-plane region , and i s negligible at polar 



latitudes. However, recent simulations by 



Tanigawa and Ikomal (120071 1 suggest a different 



picture, where accretion onto the planet mainly occurs through the vertical direction, and 
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there is no inflow along the mid-plane. In fact, iTanigawa and Ikomal (120071 ) show that there 
is significant mass outflow from the Hill sphere in the mid-plane, through the Lagrangian 
points LI and L2. This flow structure could b e a ro b ust feature of inviscid simula tions, as it 
was also observed before by iKlahr and Kleyl (120061 ). ITanigawa and Ikomal (120071 ) speculate 



that the main difference between these two results is the presence of explicit viscosity, 
as we have in the simulations presented in this work. In this case, the circumplanetary 
disk has an inflow of gas and and outflow of angular momentum. However, contrary to 
circumstellar disks, we find that the outer circumplanetary disk is yery sub-Keplerian. 



This feature of the disk has also been observed in pre vious studies ( ITanigawa and Ikoma 



2007 



Martin and Lubow 



2011 



Ayliffe and Bate 



20091 ). Unfortunately, we cannot yet 



achieve such high resolution as in their work, to properly resolve the very inner planet disk 
(l r — r p\/ r h < 0.1) where rotation might follow a Keplerian profile. This transition from 
sub-Keplerian to Keplerian rotation can have important implications for the dynamics of 
solids in circumplanetary disks. In the outer parts, large particles that orbit with Keplerian 
speed would feel a large head wind which leads to fast inwards migration. This migration 
could be significantly reduced at the location where the disk becomes Keplerian, creating 
the possibility for accumulation and growth of solids in this transitional radius. This is a 
issue that deserves further study. 

The inclusion of the curvature, global gradients, gap opening and a locally isothermal 
equation of state (all which are included in this study) could also play a role in the difference 
in flow structu re as compared to previous studies, since they break the symmetry observed 



in the work by ITanigawa and Ikomal (120071 ) . We have studied the scenario where the stellar 
disk turbulence indirectly affects the accretion of gas in the circumplanetary disk, but we 
have not assumed that the circumpl anetary disk is MRI-unstable, since the pla net disk may 



or may not be subject to the MRI ( iFujii et al 



2011 



Lubow and Martin 



20121 ). and might 



in fact have a different viscosity and accretion mechanism than the circumstellar disk. 
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Fig. 13. — Mass accretion rates by the planet for different values of a(crosses) and the 
turbulent run (diamond and dotted line). The diamond corresponds to the global average 
measure of olmhd in the magnetized disk simulation. However, this quantity varies vertically 
and decreases towards the mid-plane. 

5. SUMMARY 



For the first time, we study the accretion rate across a gap, onto Jupiter-mass planets, 
in stellar disks that have turbulence present due to the magneto-rotational instability. 
We compare our results with the planet accretion rates found in laminar viscous disks, 
and with results from previous studies, where magnetic fields and turbulence are not 
explicitly present. Contrary to all previous accretion studies, the turbulence in the disk is 
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self-consistently generated by the MRI. We also focus exclusively on the stage when the 
gap is cleared. This is a critical stage, since is a factor that determines the limiting mass 
of giant planets. We also study the structure of the outer parts of the circumplanetary 
disk that forms around the planet and give estimates that can be implemented in planet 
population synthesis models. We find that the accretion flow into the Hill sphere of the 
planet is not spherically or azimuthally symmetric, and is predominantly restricted to the 
mid-plane region of the stellar disk. Even in the turbulent case, we find no significant flow 
of mass into the Hill sphere coming from high latitudes. Accretion rates are most closely 
approximated analytically by using the reduced density in the gap region. This means that 
the gap-opening planet never reaches an accretion rate as high as the one given by the 
unperturbed density of the circumstellar disk. In a turbulent magnetized stellar disk with 
initial global stress parameter of olmhd = 1 x 10~ 3 , we find a third of the accretion rate onto 
the planet measured in a laminar viscous disk with a = 1 x 10~ 3 . The accretion rate for a 
Jupiter mass planet in an MRI-turbulent disk, after gap opening, is M 2 x 10~ 6 Mj/yr 
for a gap surface density of YlgcmT 2 . This means the growth timescale is around a few 
million years to double the mass of a Jupiter mass planet. The density profile in the outer 
parts (|r — r p \/r h > 0.3) of the circumplanetary disk can be approximated by ~ r ~ 3 / 2 and 
the rotation of the gas is highly sub-Keplerian, which has important consequences for the 
dynamics of solids in the disk around the planet and the formation of satellites. 

The authors acknowledge the computing time on the Bluegene/P supercomputer and the 
THEO cluster at the Rechenzentrum Garching (RZG) of the Max Planck Society. A. Uribe 
would like to acknowledge the support of the International Max Planck Research School 
(IMPRS) of Heidelberg. 



-27- 

REFERENCES 

Alibert, Y., Mordasini, C, and Benz, W.: 2004, A&A 417, L25 

Alibert, Y., Mordasini, C, Benz, W., and Winisdoerffer, C: 2005, A&A 434, 343 

Ayliffe, B. A. and Bate, M. R.: 2009, MNRAS 397, 657 

Ayliffe, B. A. and Bate, M. R.: 2009b, MNRAS 393, 49 

Balbus, S. A. and Hawley, J. F.: 1991, ApJ 376, 214 

Bate, M. R., Lubow, S. H., Ogilvie, G. L, and Miller, K. A.: 2003, MNRAS 341, 213 

Benz, W., Mordasini, C, Alibert, Y., and Naef, D.: 2008, Physica Scripta Volume T 
130(1), 014022 

Boss, A. P.: 1997, Science 276, 1836 

Bryden, G., Chen, X., Lin, D. N. C, Nelson, R. P., and Papaloizou, J. C. B.: 1999, ApJ 
514, 344 

D'Angelo, G., Kley, W., and Henning, T.: 2003a, in D. Deming & S. Seager (ed.), Scientific 
Frontiers in Research on Extrasolar Planets, Vol. 294 of Astronomical Society of the 
Pacific Conference Series, pp 323-326 

D'Angelo, G., Kley, W., and Henning, T.: 2003b, ApJ 586, 540 

Durisen, R. H., Boss, A. P., Mayer, L., Nelson, A. F., Quinn, T., and Rice, W. K. M.: 2007, 
Protostars and Planets V pp 607-622 

Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., and Henning, T.: 2010, A&A 515, 
A70+ 

Flock, M., Dzyurkevich, N., Klahr, H., and Mignone, A.: 2010, A&A 516, A26+ 



-28 - 

Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., and Henning, T.: 2011, ApJ 735, 122 

Fujii, Y. L, Okuzumi, S., and Inutsuka, S.-i.: 2011, ApJ 743, 53 

Gardiner, T. A. and Stone, J. M.: 2005, Journal of Computational Physics 205, 509 

Hubickyj, O., Bodenheimer, P., and Lissauer, J. J.: 2005, Icarus 179, 415 

Ida, S. and Lin, D. N. C: 2008, in Y.-S. Sun, S. Ferraz-Mello, & J.-L. Zhou (ed.), IAU 
Symposium, Vol. 249 of IAU Symposium, pp 223-232 

Ikoma, M., Nakazawa, K., and Emori, H.: 2000, ApJ 537, 1013 

Klahr, H. and Kley, W.: 2006, A&A 445, 747 

Kley, W., D'Angelo, G., and Henning, T.: 2001, ApJ 547, 457 

Lubow, S. H. and Martin, R. G.: 2012, ApJ 749, L37 

Lubow, S. H., Seibert, M., and Artymowicz, P.: 1999, ApJ 526, 1001 

Martin, R. G. and Lubow, S. H.: 2011, MNRAS 413, 1447 

Mayer, L., Quinn, T., Wadsley, J., and Stadel, J.: 2002, Science 298, 1756 

Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C, and Ferrari, 
A.: 2007, ApJS 170, 228 

Mizuno, H.: 1980, Progress of Theoretical Physics 64, 544 

Mordasini, C, Alibert, Y., and Benz, W.: 2009a, A&A 501, 1139 

Mordasini, C, Alibert, Y., Benz, W., and Naef, D.: 2009b, A&A 501, 1161 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., and Greenzweig, 
Y.: 1996, Icarus 124, 62 

Rafikov, R. R.: 2005, ApJ 621, L69 



-29 - 

Shakura, N. I. and Sunyaev, R. A.: 1973, A&A 24, 337 

Stevenson, D. J.: 1982, Planet. Space Sci. 30, 755 

Tanigawa, T. and Ikoma, M.: 2007, ApJ 667, 557 

Tanigawa, T., Ohtsuki, K., and Machida, M. N.: 2012, ApJ 747, 47 

Tanigawa, T. and Watanabe, S.-i.: 2002, ApJ 580, 506 

Turner, N. J., Carballido, A., and Sano, T.: 2010, ApJ 708, 188 

Uribe, A. L., Klahr, H., Flock, M., and Henning, T.: 2011, ApJ 736, 85 

Ward, W. R. and Canup, R. M.: 2010, AJ 140, 1168 

Winters, W. F., Balbus, S. A., and Hawley, J. F.: 2003, ApJ 589, 543 



This manuscript was prepared with the A AS IAT^X macros v5.2. 



